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Abstract 

, _ _ , This paper addresses the problem of correlation estimation in sets of compressed images. We consider a framework where 

^-H . images are represented under the form of linear measurements due to low complexity sensing or security requirements. We assume 

^-H ■ that the images are correlated through the displacement of visual objects due to motion or viewpoint change and the correlation 

C^ ' is effectively represented by optical flow or motion field models. The correlation is estimated in the compressed domain by 

CS| , jointly processing the linear measurements. We first show that the correlated images can be efficiently related using a linear 

- X . operator. Using this linear relationship we then describe the dependencies between images in the compressed domain. We further 

^ ■ cast a regularized optimization problem where the correlation is estimated in order to satisfy both data consistency and motion 

Q smoothness objectives with a Graph Cut algorithm. We analyze in detail the correlation estimation performance and quantify the 

\ penalty due to image compression. Extensive experiments in stereo and video imaging applications show that our novel solution 

QQ . Stays competitive with methods that implement complex image reconstruction steps prior to correlation estimation. We finally use 

^_^ ■ the estimated correlation in a novel joint image reconstruction scheme that is based on an optimization problem with sparsity 
priors on the reconstructed images. Additional experiments show that our correlation estimation algorithm leads to an effective 

1—^ , reconstruction of pairs of images in distributed image coding schemes that outperform independent reconstruction algorithms by 

^ . 2 to 4 dB. 

^^ I Index Terms 

(^ _ Linear measurements, correlation estimation, distributed image compression, joint reconstruction 

CS| ■ L Introduction 

J^ . T^ recent years, the increasing popularity of vision sensor networks has led to the generation of huge volume of visual 
\^ ■ i information. This creates the need for effective information processing systems that are able to efficiently compress, analyze 
\^ and store highly redundant information streams captured by multiple devices. Distributed processing solutions become highly 
"^ attractive in such a context, as they permit to reduce the communication and computational power requirements in the sensors. 
f^^ The visual information is typically compressed and transmitted independently from each sensor node to a common decoder that 
^^ ' jointly processes the correlated information streams. The inter-sensor communication needs are relaxed and the computational 
^"^ burden is shifted to the decoder. The estimation of the image correlation at decoder becomes however crucial in such distributed 
settings for image reconstruction or analysis tasks. 

In this paper, we consider the problem of correlation estimation in a framework where multiple sensors transmit compressed 
images that have been obtained by a small number of linear projections of the original images, as illustrated in Fig. [T] Such 
linear projections typically represent simple measurements in low complexity sensing systems El, 0. We propose a novel 
C^ ' solution for correlation estimation at the joint decoder where the analysis is performed directly in the compressed domain in 
order to avoid expensive image reconstruction tasks. This is especially useful for the analysis applications that do not target 
image reconstruction. We assume that the correlation between images corresponds to camera or object motion; this can be 
efficiently represented by optical flow or motion field models. We show that such a correlation model can be described by 
a linear operator and we further analyze in detail the effect of such operator in the compressed domain. Later, we cast the 
correlation estimation as a regularized energy minimization problem with constraints on data consistency as well as consistency 
of the motion field. In particular, we regularize the correlation model such that the motion values in neighboring pixels are 
similar except at image discontinuities. Such an optimization problem can be solved by Graph Cuts algorithms. 

We analyze in details the performance of our novel correlation estimation framework. In particular, we study the penalty 
in the correlation estimation that is due to working in the compressed domain as opposed to the original image domain 
as in traditional correlation estimation problems. We show that the penalty decreases when the number of measurements 
increases and that our algorithm tends to the optimal correlation estimate at high measurement rate. Extensive simulations in 
distributed stereo and video imaging applications confirm that the proposed solution provides effective estimates of the relative 
motion between images and even competes with solutions that implement expensive image reconstruction prior to correlation 
estimation. We finally study the performance of a novel joint reconstruction algorithm that uses our correlation estimates for 
decoding pairs of images. The joint reconstruction is cast as an optimization problem where the reconstructed images have 
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Fig. 1. Schematic representation of the proposed scheme. The images /i and I2 are correlated through displacement of scene objects due to viewpoint 
change or motion of scene objects. The correlation model is estimated directly in the compressed domain without any intermediate image reconstruction. The 
correlation information can be used for optional joint reconstruction. 

to satisfy sparsity priors as well as consistency with both the measurements and the correlation estimates. We solve this joint 
reconstruction problem by effective proximal splitting methods and show that accurate correlation estimation in distributed 
image representation permits to outperform independent decoding solutions in terms of image quality. 

The rest of the paper is organized as follows. Section [III briefly overviews the related work. In Section |III1 we describe 
the proposed framework and show how the correlation estimation problem carries out to the compressed domain. Section |IV] 
describes the proposed correlation estimation algorithm and its performance are analyzed in details on Section |Vl In Section 
rvJ we draw some concluding remarks. 

II. Related WORK 

This section describes the literature related to the framework proposed in this paper. We first present sensing solutions based 
on linear measurements. We then review the correlation estimation algorithms in distributed image representation systems and 
finally we discuss the most relevant works about joint reconstruction of correlated images. 

In recent years, signal acquisition based on random projections received a significant attention in many applications like 
medical imaging, compressive imaging and even sensor networks. Donoho [2| and Candes et ah |3| show that the small number 
of linear measurements contain enough information to reconstruct a sparse or a compressible signal. In particular they show that 
if a signal has a sparse representation in one basis then it can be recovered from a small number of linear measurements taken 
on another (random) basis that is incoherent with the first one. Essentially, if the signal is iiT-sparse (i.e., if the signal contains 
K significant components), then one need approximately cK linear measurements (typically c = 3 or 4) to reconstruct the 
signal with high probability [4] . Such results open the door to novel low complexity sensing solutions where the computational 
complexity for signal reconstruction or analysis is pushed to the decoder. These ideas have been applied to image acquisition 
0, in, Q and later extended to video sequences lUI, (91, ifTOl . ifTTl . The effect of measurement quantization and the lossy 
compression of linear measurements has been studied in |[T2l . 

One of the key characteristics in imaging applications resides in the high correlation between multi-view images or successive 
images in a video sequence. The correlation could be exploited for effective reconstruction of image sets or for joint analysis 
tasks in distributed systems. Duarte et ah |[T3]| , llT4t have proposed different correlation models for the distributed compression 
of correlated signals from linear measurements. In particular, they introduce three joint sparsity models (JSM) in order to 
exploit the inter-signal correlation for the joint reconstruction. They are respectively described by (i) JSM-1, where the signals 
share a common sparse support plus a sparse innovation part specific to each signal, (ii) JSM-2, where the signals share 
a common sparse support with different coefficients, and (iii) JSM- 3 with a non- sparse common signal with an individual 
sparse innovation in each signal. These correlation models permit an effective joint reconstruction with a small number of 
measurements compared to independent reconstruction. These simple joint sparsity models are however not ideal for multi-view 
images or video sequences, as the correlation model in such scenarios is usually given in the form of disparity or motion vectors 
respectively. The authors in ifTSll . O, ifTTl have proposed a distributed joint reconstruction scheme for video sequences based 
on linear measurements. These schemes split the video sequences into key frames and compressed sensing (CS) frames. The 
random measurements are computed independently for each compressed sensing frame and are transmitted to the joint decoder. 
The key frames are intra coded and the joint decoder builds the side information from the intra coded key frames. The generated 
side information is then used to decode the CS frame by solving an optimization problem which assumes that the prediction 
error for the CS frame is sparse in an orthonormal basis |[T6l or block-based adaptive dictionary ifTTIl . Similar ideas have been 
used by Trocan et al. ifTSl . |[T9l for distributed multi-view compression where all images are however given in the form of 
linear measurements. The joint decoder first reconstructs all the views by solving a regularized optimization problem. Then, 
the independently reconstructed views are used to estimate the underlying correlation model in the form of disparity image. 
The disparity image is used to jointly reconstruct all the views in a multistage refinement framework in which each refinement 
stage reconstructs a view by assuming that the prediction error is sparse in a dual tree wavelet basis. In our previous work ll20l . 
we have proposed a methodology to estimate the correlation between pairs of frames where one image serves as a reference 
image. Unfortunately, reconstructing the reference image in a compressed measurements framework typically requires methods 



based on solving l2-h optimization problems that are highly complex. The works in the literature typically use reference 
frames that are encoded as intra- frames or reconstructed with complex optimization tools prior to correlation estimation. In 
this paper, we rather propose to avoid the explicit reconstruction of the images and directly estimate the correlation in the 
compressed domain. In general, it permits to reduce the computational complexity at decoder, especially in applications where 
image reconstruction is not necessary. 

Finally, other works have recently addressed the problem of joint reconstruction of correlated images given in compressed 
form. For example. Park et al I2JL1 have proposed an image registration and joint reconstruction algorithm for multi-view 
images based on manifold lifting. The underlying correlation between views is exploited by constructing an image appearance 
manifold where the images represent sample points on the manifold; these points are controlled by a few camera parameters 
(e.g., rotation, translation etc.). By knowing the initial camera positions the images are jointly reconstructed based on an I2-I1 
optimization framework. Then, the reconstructed scene is used to refine the camera parameters and thus both the camera 
positions and the scene are jointly estimated using alternating minimization techniques. In another framework 12211 . the authors 
have proposed a joint reconstruction scheme based on a regularized optimization framework. The two regularization terms 
encourage sparse priors of multi-view images and their difference images. However, the correlation between images is not 
efficiently exploited as the correlation model in multi-view image settings is usually given in the form of disparity image and not 
as a sparsity prior of the signal differences. The joint reconstruction scheme proposed in this paper rather builds the correlation 
model in the form of disparity or motion field and thus facilitates an efficient joint image representation. Furthermore, the 
correlation model is built directly from the linear measurements in the compressed domain and thus avoids the computational 
complexity of reconstructing the reference images. The proposed scheme is based on low complexity linear measurements 
and it provides an interesting flexible solution for distributed processing in vision sensors, targeting applications like object 
detection, distributed rendering or distributed joint signal reconstruction. 



III. Distributed Representation of Correlated Images 
A. Framework 

We consider a framework where the images represent a scene at different time instants or from different viewpoints. For 
the sake of clarity, we consider a pair of images /i and I2 (with resolution N=Ni x N2) but the framework extends to larger 
number of images. These images are represented by linear measurements that correspond to the projection of the image pixel 
values on a set of coding vectors. Typically, the coding vectors can be constructed from Gaussian or Bernoulli distributions ^ 
or with a block structure (Tl, 1231 for easier handling and fast sampling of large images. The measurements are transmitted to 
a joint decoder that can estimate the correlation between the compressed images and possibly perform a joint reconstruction 
of the image set. The framework is illustrated in Fig. [T] 

In more details, the sensors process images row by row. Let Ii^k and /2,/e represent the k^^ row of the images /i and I2 
respectively, and Yi^k and F2,fe represent the linear measurements computed from Ii^k and /2,fe using the measurement matrices 
(l)i and (/)2 respectively. The measurements Yi^k and F2,fe are computed as 

(1) 

,2,...,iVi, 

where (.)^ denotes the transpose operator. It should be noted that (f)^ and 02 are of dimensions M x N2, where M « N2 
is the number of measurements computed for each row in the image. From Eq. ([T]) it is easy to check that the measurements 

^1 = [^1,1, ^1,2, • • . yi,Ni] and I2 = [^2,1, ^2,2, . . • y2,Ni] can be computed as 
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where $, is the measurement matrix used to sample the i*'^ image Vi = {1, 2}. It is represented as 
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where K = MNi, N = N1N2 and K/N represents the measurement rate. 
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Fig. 2. Illustration of the relation between the matrices A and B. On the left the images /i, /2 G M^ are related using the matrix A where A^ = N1N2. In 
the compressed domain, the measurement vectors Yi , 12 ^ '^^ are related using the matrix B where K = MNi . The matrices A and B can be related by 
B ^ $2^<^{ where <l>i's are the sensing matrices. 

B. Correlation Model 

In the above settings, the correlation between images can be mainly explained by the relative displacement of objects in the 
scene. This can be modeled effectively by the optical flow that determines the amount of displacement of objects or pixels in 
different images. We show now, how such a correlation model can be described by a linear operator. Let m^ and m^ represent 
the horizontal and vertical motion components. As the visual objects in the images /i and I2 are displaced, it is likely that the 
pixel at position z = (/c, /) in one image moves to z' = {k -\- m^(/c, /), / + m^(/c, I)) in the second image. Thus, the images 
/i and I2 can be simply related by a linear operator T that changes the coordinate system from (/c,/) in the first image to 
{k + m^(/c, /), / + m'^(/c, /)) in the second image, i.e.. 



h = T{h} 
For mathematical convenience we use an equivalent representation of Eq. (|4]) in the form of matrix multiplication: 
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2,/e 
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^l,Ai 



V^ = l,2,...,7Vi, 



(4) 



(5) 



where A^ is a matrix of dimensions N2 x N1N2 whose entries are determined by the horizontal and vertical components of 
the motion field in the k^^ row of pixels, i.e., m^(/c, .) and m^(/c, .). The elements of the matrix A^ are given by 



A^(/,/ + /3i+/32iV2) 



_ / 1 if m^(A;,/) = /3i,m^(/c,0 =/32 
otherwise 



(6) 



If I -\- Pi -\- I32N2 > N1N2 (e.g., at image boundaries), we set I -\- Pi -\- I32N2 = N1N2 so that the dimensions of the matrix 
A^ stays N2 x N1N2. It is easy to check that the matrix A^ formed using Eq. (O contains only one '1' in each row; this 
implies h^kil) = h,k+(3i{l -^ h) if A^(/, / + /3i + /32N2) = 1. The action of the matrix A^ shifts the pixels in /i according 
to the motion given as m^(/c, .) and m^(/c, .) and forms an estimate of the image /2,/c- It should be noted that the matrix A^ 
is completely determined by the k^^ row of motion vectors m^(/c, .) and m^(/c, .). 

The relation given in Eq. (O can be extended to all rows of the image I2. The images /i and I2 are finally related by a 
linear operator A such that I2 = A Ii which can be written as 



(7) 



r T^ 1 
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h A 

This relation is illustrated on the lefthand side of Fig. |2l 



h 



C. Correlation in the measurement domain 

We now extend the above correlation model in the compressed domain. Without loss of generality, we first assume that the 
measurements Yi and I2 can be related by a linear transformation B, i.e., 



(8) 
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where B^, \/k = 1, 2, . . . , TVi is a matrix with dimensions M x MNi, i.e., the measurements l2,/c can be related to Yi as 

1^2,fc = 5^yi, Vfe = l,2,...,iVi. (9) 

Any two vectors Yi , F2 ^ M^^^ can be related by a linear transformation 5 as long as Fi 7^ 0, which is the case in our 
framework. However, this linear transformation B a priori does not have any special form. We are interested in understanding 
the relation between this matrix and the matrix A that shifts the pixels between images h and I2. Pre-multiplying Eq. (|7]) by 
<l>2 on both sides, one can write 

^2 = ^2/2 = ^2^/1. (10) 

In addition, by replacing Yi = ^ili in Eq. ([8]), one obtains 

Y2 = BYi = B^ih. (11) 

From Eqs. JTOl) and (fTTt the relation between B and A can finally be given as 

B^i = <I>2A (12) 

This forms an over-determined system of linear equations, as the number of unknown in matrix B is smaller than the number 
of equations in $2^. In this case, the optimal matrix B that minimizes \\B^i — <l>2^||2 is given by 



B 



^2^^l^ 



(13) 



where f denotes the pseudo-inverse operator. As the rows in <l>i are generally orthonormal, the pseudo-inverse can be computed 

using the 

becomes 



using the transpose operator, i.e., <l>i^ = ^1^. Substituting B = ^2^^!^ in Eq. ^ the relation between the measurements 
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(14) 
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Eq. ([T5]) comes from the fact that measurements are computed across rows of pixels in our framework. The relationship between 
the matrices A and B is illustrated in Fig. [21 where the matrices A and B are used to relate the points in the original and 
compressed domains, respectively. In the next section, we propose an algorithm for estimating the correlation model directly 
from the linear measurements Yi and Y2. 

IV. Correlation estimation from linear measurements 

A. Regularized energy minimization problem 

We propose in this section a method for estimating the correlation between images from the compressed measurements 
without any explicit image reconstruction step. The objective is to compute a flow or motion field that represents the motion 
between images /i and I2. We denote this flow field as A^ = (m^,m^), where m^ and m^ are horizontal and vertical 
components of the motion, respectively. The problem then consists in finding the value of the flow field M at each pixel 
position z = (/c, /) such that the estimated correlation is consistent with the measurement vectors Yi and Y2. At the same 
time, the motion field has to be piecewise smooth in order to model consistent motion of visual objects. We propose to cast 
the correlation estimation as a regularized energy minimization problem where the energy E{M) is composed of a data term 
Ed{M) and a smoothness term Es{M). The optimal mapping A^* is obtained by minimizing the energy function E{M) as 

A^* = arg min E{M) = arg mm[Ed{M) + XEs{M)], (16) 

M M 

where A balances the importance of the data and smoothness terms. 

We now discuss in more details the components of the energy function. The smoothness term measures the penalty of 
assigning different motion values to the adjacent pixels. We write it as 

Es{M)= YI VMiz^z), (17) 



where z, il are neighbour pixels in the 4-pixel neighbourhood denoted by M . The term V}v4(z, z') is given as 

ia^(z,zO = mm(|m^(z)-m^(zO| + |m^(z)-m^(zO|,r), (18) 



where r sets an upper level on the smoothness penalty that helps preserving the discontinuities 

Next, the data function measures the consistency of a particular motion value for pixel z with the vectors Y\ and Y2- 
Classically, the accuracy of the motion values is evaluated with the original images ll25l and the data cost is typically given as 

Ed{M) = Y,Y.^M{k,l), (19) 

where AM{k^ I) = \\h,k{l) — A,fe+m^(fe,o(^ + m^(/c, /))||2 represents the error of matching the pixel at position (/c, /) in the 
second image with a pixel in the first image that is selected according to the motion information. As discussed in Section 
IIII-BI the effect of motion between images can be captured by a linear operator A that is a composition of sub-matrices A^ . 
We can therefore rewrite Eq. ([T9l ) as 

Ed{M) = Y^\\ll,-AX,h\\l (20) 

where the sub-matrix A^ depends on the motion field M according to Eq. (|6]). In the rest of the development, we drop the 
index A4 as the dependency on the motion field is clear from the context. 

In our framework however we do not have access to the original images, but only to the measurement vectors Yi and Y2. 
We thus approximate the data cost Ed{M) by Ed{M) that is computed directly from the measurement vectors. It can be 
written as 

EdiM) = Ell^2':.-^'^i|l2 (21) 

k=l 

Ni 

^ ^\\Y2,k-B^Y,\\l (22) 

k=i 



^\\Y2,k-^'2A'^i^Y,\\l (23) 



= Ed{M). (24) 

Note that the data cost approximation due to working in the compressed domain comes from the relation between the matrices 
A and B that is given as 5 ^ ^2^^! • We study in details the effect of this approximation in the next section. 
We can finally rewrite the regularized energy objective function for the correlation estimation problem. It reads as 

E{M) = Y,\\Y2,k-^^2Aii^i^Y,\\l^X ^ Vm{z^z). (25) 

/c=l z,z'GA/' 

This cost function is used in the optimization problem of Eq. ([T6l) , which becomes a non-convex problem. The search space 
is discrete and usually constrained by limits on each of the motion values which typically define a motion search window. 
The solution to this problem can be determined with strong optimization techniques based on Graph Cuts ll24l . l|26l or Belief 
propagation fTf\. A comprehensive overview of various energy minimization techniques is summarized in ll28ll . In this paper, 
we use an optimization algorithm based on a-expansion mode in Graph Cuts whose complexity is bounded by a low order 
polynomial ll26l . 

Finally, it should be noted that the correlation estimation can also be performed by block of pixels. In this case, each block 
of pixels is assumed to move in a coherent way, and the objective of the correlation estimation problem is to compute one 
motion vector per block. The data cost function can then be modified in a straightforward way by imposing the same motion 
vector for all the pixels in a block. The smoothness function is also modified in this case such that it penalizes the difference 
between the motion values of adjacent blocks rather than neighboring pixels. The optimization problem keeps the same form 
in block-based motion estimation but the search space is dramatically reduced as the number of motion vectors is smaller. 

B. Compressed domain penalty 

We now discuss the penalty of estimating the correlation from measurements instead of original images. When the correlation 
M is given, this penalty corresponds to the difference between the values of the regularized energy function of Eq. ([T6l) that is 
evaluated from original images or respectively measurements. Recall that the smoothness cost function Es{M) depends only 
on the correlation (see Eq. ([TT])). Therefore, the estimation penalty is identical to the error between the data cost functions 



Ed{M) and Ed{M) that are computed in the original and compressed domains respectively. We first show that the penalty is 
bounded. Then, we show that the penalty decreases monotonically when the number of measurements increases. 

Proposition 1: The penalty of estimating the correlation from measurements is bounded. In particular we have | (1 — 5)'^E(i{M.)- 
Ci\< Ed{M) < {I ^5fEd{M)^Cu. ^\i^rt 5 >{), Ci = T]^ ^2{l- 5)ar]^ Cu = T]^ ^2{l^ 5)a^^ 

ri = Ek=l ^maM') T!;t7-^, Wkv - hA2 with h = ^f 1^. 

Proof: Let us assume that M or equivalently A is given. Then, the points V = {/J/^, A^Ii : k = 1, 2, • • • , A^i} forms 
a finite set. According to the Johnson-Lindenstrauss (JL) lemma, the distances between points in V are preserved in the 
measurement domain R^ when M = 0{S~'^ log\V\) measurements are computed with a measurement matrix ^2 lE9l , ll30l , 
where \V\ denotes the number of points in V. Mathematically, the JL-embedding is given as 

(1 - 5)\\ll, - A'^hh < U'llk - 02'^'/ill2 < (1 + S)\\ll, - A'^hh, (26) 

for a positive constant S. It should be noted that, when the measurement matrix 02 satisfies Eq. (|26l) , then with high probability 
it satisfies the restricted isometry property (RIP). For more details related to the connection between the JL-lemma and the RIP 
we refer the reader to l|29l , ll30l . Eq. (l26l) holds with high probability not only for measurement matrices constructed using 
Gaussian and Bernoulli distributions but also for structured measurement matrices constructed using orthonormal bases, e.g., 
DCT, EFT [30|. In our experiments we construct measurement matrices using structured EET. 
Eor a given row index k, the term F2,/c — (/>2^^<^f Yi in Eq. (|23]) we can write as 

= ^lllu-^lA^^^c^lA^h-c^lA^^l^rh 

= ^2llk-^2A^h^Eh, (27) 

where E = 02^^ - 02^^*^r*^i- The term ||F2,/c - (/>2^^^f ^ilb can be upper bounded as 



\Y2,k-^'2A'^lY,\\2 = U'2llk-^2A'h^Eh\ 



2 



< {l + 5)\\ll,-A'^hh + \\Eh\\2, (28) 

where the last inequality is derived from Eq. (l26l) . Similarly the term \\Y2^k — 02^^<^fyi||2 can be lower bounded as 

\\Y2,k - 0t^'*r^i II2 = U^2llk - <PlA^h + Eh II2 

= Ulllu-<PlA''h-{-Eh)h 

> lU'llk-^A'^hh-WEhhl (29) 

> 1(1- ,5)11/2^, -A'=/i||2-||^/l||2|, (30) 

where Eq. ^ follows from \\x - y\\2 > \\\x\\2 - ||y||2|, and Eq. ^ is derived from Eq. (Egji. The term ||£;/i||2 in Eqs. ^ 
and (l30l ) can also be bounded as 

\\Eh\\2 = ||(/.^A'=Ji-(/.^A'=$f$iJi||2 

= \\4A\^J^ih-h)\\2 

< \\A\h-h)\\2 (31) 

k-\-Wy 

< \\A% J2 Il^i.p-^i.pll2 

P—k—Wy 

k-\-Wy 

= (TmaxiA'') Y^ \\h,p - h,ph = r]k, (32) 

P^k—Wy 

where Eq. (ISTt follows from ||<^x||2 < ||^||2, as $ is a non expanding operator f29l and /i = <l>f $1/1 is the pre-image 
of /i. (Jmax{A^) in Eq. (l32l) denotes the largest singular value of A^ . The summation in Eq. (l32l) is carried out from rows 
k — Wyiok^Wydi^ the search window is usually bounded, where Wy is the admissible search size along the vertical direction. 
Combining Eq. (|28l) , Eq. (l3Ql) and Eq. (l32l) , and by taking squares we get for each row of pixels 

{\{i-5)\\ii^-A^i42-vk\f<\\y2M-<PlA'^^TyAl<{{^ + milk-A^hh+vkf. (33) 

Adding the second and third inequality terms of Eq. (|33]) for all rows /c = 1, 2, • • • , TVi results in 

A^i 

MM) = Y.\\Y2^k-^lA^^lY\\l 

k=l 



where 



and 



\k=l ) fe=l fe=l 



= (1 + SfEd{M) + r?2 + 2(1 + b)r]a 

= (l + J)2^rf(7W) + C„, (34) 

C„ = 7?2 + 2(1 + ,5)7?a, (35) 



iVi iVi k^Wy 

r] = ^r]k = ^(y^ax{A^) X! Il^i,p-^i,pll2- (36) 

/c=l /c=l p—k—Wy 

In a similar way, from the first and second inequahty terms of Eq. ([33]) we get 



k=i 



Ni 



> 






■ ATi \ 2 ATi ATi 



where 



and T] is given in Eq. ([36b . 



>|(l-5)2^,(A1)-K],?J -2(l-5)^,?fc^||/2^,-A'=/i||2| 

77 a 

= |(l-(5)2E,(M)-(7/2 + 2(l-(5)7?a)| 

= \{l-6fEd{M)-Cil (37) 

a=7^2 + 2(l-(5)7?a, (38) 



Proposition 2: The penalty of estimating the correlation from measurements monotonically decreases when the measurement 
rate K/N increases. It further becomes negligible at high measurement rate. 

Proof: In Proposition [T] we have shown that the difference between the data cost functions estimated from compressed 
measurements Ed{M) and images Ed{M) is lower and upper bounded by errors Ci and Cn, respectively given in Eq. ([38]) and 
Eq. ([35]) . The error r] = Xlfe=i ^^ ^ Xl/c=i ll^i,^ ~ ^i,fe||2 (see Eq. ([36]) ) decreases with increasing measurement rate because 
{(fy^ (j)^ becomes an orthogonal projection operator, and h = <l>f <l>i/i becomes arbitrarily close to h when the number 
of measurements increases. Therefore, the errors Ci and Cu decrease as the measurement rate increases and when sufficient 
number of measurements are taken the errors Ci and Cu become negligible, i.e., E(ji{M) ~ Ed{M). ■ 

Due to the error between the cost functions Ed{M) and Ed{M), the solution M estimated from the linear measurements is 
not accurate especially at low measurement rates. The solution of the correlation estimation problem in the compressed domain 
might thus be quite far from the actual correlation between images. However, as the number of measurements increases the 
approximation in the compressed domain becomes more accurate and the solution of the correlation estimation tends to the 
actual correlation between original images. 



V. Experimental results 

A. Setup 

We analyze the performance of the correlation estimation in both stereo and video imaging applications. The random 
projections are computed using a scrambled Fourier measurement matrix where the scrambled operator is a diagonal matrix 
with entries ±1 taken from an i.i.d. Bernoulli random variable with equal probability |30|. In all our experiments, we sample 
both images using a same measurement rate. The correlation is estimated without prior reconstruction of images by minimizing 
the objective function in Eq. ([25l) . 

For the stereo imaging case, we evaluate the disparity estimation performance in two natural image sets namely Tsukuba 
and Venu^ 1311 . These datasets have been captured by a camera array where the different viewpoints correspond to translating 
the camera along one of the image coordinate axis. In such a scenario the motion of objects due to the viewpoint change 
is restricted to the horizontal direction with no motion along the vertical direction. The disparity estimation is thus a one- 
dimensional search problem, and the smoothness and data cost functions are modified accordingly by assuming that m^ = 0. 
The size of the search windows used in our experiments are 16 pixels for Tsukuba and 20 pixels for Venus. In our experiments 
we estimate disparity in both dense (per pixel) and block settings, where the block size is fixed to 4 x 4 pixels. 

In the video scenario, we analyze the motion estimation accuracy in two synthetic scenes, namely Yo Semite and Grove, and 
one natural scene Mequoil^. The Grove and Mequon datasets are resampled to a resolution of 160 x 120 pixels using bilinear 
filters. The size of search windows is of ±3 pixels in both horizontal and vertical directions. For the sake of simplicity, we 
estimate a motion field for blocks of pixels with size of 4 x 4 pixels and not motion vectors for each pixel. 

The accuracy of the correlation estimation is evaluated by comparing to groundtruth information and to the correlation 
estimated from the reconstructed images. We propose another representation of the accuracy of the correlation by discussing 
the quality of the second view I2 that is reconstructed by prediction of the first frame according to estimated correlation. We 
then analyze the influence of the sampling matrix and the effect of measurement quantization on the correlation estimation 
performance. We finally show the importance of accurate correlation estimation in a novel joint reconstruction algorithm where 
images are reconstructed from measurements while satisfying sparsity constraints as well as consistency with the correlation 
information. Note that in practice, the groundtruth correlation model and the original images are not available a priori to estimate 
an optimal regularization parameter A. In such cases, the regularization parameter A can be estimated based on learning from 
a set of training images or using the automated method proposed in f32|. In our experiments, we however select the parameter 
A based on trial and error experiments. 

B. Disparity estimation performance 

We first illustrate disparity maps for the Venus dataset where the compressed data have been obtained with a different 
measurement matrix for each image, i.e., <l>i 7^ <l>2. Fig. Ob) and Fig. [Sfd) show the disparity map from a measurement rate 
0.2 and 0.7 respectively. Fig. Oc) and Fig. Oe) represent the corresponding disparity errors. Comparing the results with the 
groundtruth given in Fig.Oa) we see that at low measurement rate (corresponding to 0.2) we estimate a coarse version of the 
disparity map. Quantitatively the disparity error with respect to the groundtruth is found out to be 41% when measured as the 
percentage of pixels where the absolute error is greater than one ISTIl as shown in Fig. He). At higher rate, the disparity map 
is more accurate and the disparity error drops below 11%. 




Fig. 3. Comparison of the estimated disparity image with respect to groundtruth information at measurement rates 0.2 and 0.7 in the Venus dataset. (a) 
Groundtruth disparity image M^; (b) computed dense disparity image m^ at measurement rate 0.2; (c) disparity error at rate 0.2. The pixels with absolute 
error greater than one is marked in white. The percentage of white pixels is 41%. (d) Computed dense disparity image m^ at measurement rate 0.7; (e) 
disparity error at rate 0.7. The percentage of white pixels is 10.7%. 

We then show the quality of the reconstruction of the second image that is predicted from the first image using the correlation 
estimate. When a coarse disparity map m^ (i.e., estimated at low measurement rate) is used for image prediction the resulting 
predicted image I2 is closer to I2 than h (see Fig. Ua) and Fig. [Jb) respectively). We observe that the mean square error 



^Available in 
^Available in 



http://vision.middlebury.edu/stereo/data/scenes2001/| 
htt p ://vision . middlebury. edu/flow | 
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Fig. 4. Evaluating the accuracy of disparity image in Fig. (Sjb) and Fig. (Sjd) in terms of image prediction quality for the Venus dataset. The disparity images 
in Fig.jSjb) and Fig.[3jd) are used to predict the image I2 at measurement rates 0.2 and 0.7 respectively, (a) Inverse prediction 1 — |/2 — /2I at a measurement 
rate of 0.2; (b) inverse prediction 1 — |/2 — /i| at a measurement rate of 0.2; (c) inverse prediction 1 — |/2 — ^2! at a measurement rate of 0.7; (d) inverse 
prediction 1 — |/2 — ^i| at a measurement rate of 0.7. The error is inverted, so that the white pixels correspond to no error. 
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Fig. 5. Comparison of the quality of predicted image /2 with respect to I2 and /i with and without regularization, i.e., A / and A = in Eq. 125} 
respectively in the Tsukuba dataset. The image prediction is carried out using dense disparity image. 



(MSE) between the predicted image I2 and I2 is smaller than the error between I2 and /i, which confirms the benefit of 
the disparity estimate in the prediction. When the measurement rate increases the quality of the disparity map improves and 
the quality of the predicted image I2 also improves substantially as it can be observed in Fig. (He) and Fig. Ud) where the 
measurement rate is set to 0.7. 

The next experiments illustrate the benefit of regularization in the disparity estimation problem. Fig. \5\ plots the quality of 
the predicted image I2 with and without smoothness cost (i.e., A 7^ and A = respectively in Eq. ([25l) respectively) for the 
Tsukuba dataset. It is clear that the quality of I2 is improved by enabling the regularization term in our optimization framework. 
Similar experimental finding is observed for the Venus dataset. From Fig. [5] we further observe that the quality of the predicted 
image I2 at a measurement rate 0.05 is 22.2 dB (the corresponding disparity error is 39%), which is approximately 3.5 dB 
away from the saturation point or from global minima solution due to influence of the penalty terms Ci and Cu discussed in 
Section IIV-BI As the measurement rate increases, the influence of the terms Ci and Cu decreases. As a result the quality of 
predicted image I2 increases with the measurement rate and saturates above rates > 0.5. In other words, our scheme gives 
optimal disparity solution at high measurement rate. We also carry out experiments using the same measurement matrix for 
both images, i.e., <l>i = ^2- Fig.Oa) and Fig. Ob) compare the PSNR quality of the predicted image I2 and the disparity error 
(DE), respectively with the results obtained with different measurement matrices. It is clear that the prediction image quality 
and the disparity accuracy improve when different measurement matrices are used as this brings more information from both 
images to solve the correspondence problem. 

We finally compare our disparity estimation results to a scheme that first reconstructs the images before estimating the 
disparity map. The images are reconstructed independently from the corresponding measurements by solving a convex opti- 
mization problem. We denote this methodology as disparity from reconstructed images (DFR). We have tried out two different 
reconstruction methodologies: (1) DFR-sparsity that consists in minimizing the li norm of the sparse coefficients assuming that 
the image is sparse in a particular orthonormal basis (e.g., a wavelet basis); this problem is solved using GPSR 1331 : (2) DFR- 
TV that minimizes the TV norm of the reconstructed image; this problem is solved using BPDQ toolbox 1341 . The disparity 
map is then estimated using a-expansion mode in Graph Cuts applied on the reconstructed images. Fig. [6l shows the comparison 
of the proposed scheme with the DFR-sparsity and DFR-TV schemes for the Tsukuba dataset. From Fig. Oa) and Fig. Ob) 
we observe that the performance of our low complexity correlation solution competes with DFR-sparsity scheme. Especially at 
rates smaller than 0.1, our scheme performs better than the DFR-sparsity scheme as the poor image reconstruction quality in the 
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Fig. 6. Comparison of the proposed scheme with the DFR schemes in the Tsukuba dataset: (a) comparison in terms of image prediction quahty; (b) comparison 
in terms of disparity error. The performance of the proposed scheme is evaluated using both the same and different sets of measurement matrices. 
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Fig. 7. Comparison of the proposed scheme over disparity from reconstructed image (DFR) schemes for Tsukuba dataset. The disparity is estimated per block 
with a block size 4x4. The performance of the proposed scheme is evaluated using both the same and different sets of measurement matrices. 



DFR-sparsity scheme leads to a bad estimation of disparity map. On the other hand, DFR-TV scheme significantly outperforms 
our scheme at lower rates due to good reconstruction image quality. However, our scheme estimates the optimal correlation 
for rates above 0.5 and thus performs similar to the DFR-TV scheme at high measurement rate but with a complexity that is 
dramatically smaller as it avoids image reconstruction. In particular, in our experiments we have observed that the running time 
of Graph Cuts algorithms that estimate the correlation information from linear measurements is approximately same as the 
one that estimates the correlation information from reconstructed images. The complexity of our correlation estimation scheme 
stays reasonable due to the efficiency of Graph Cuts algorithms whose complexity is bounded by a low order polynomial 
ll26l . Comparing to the DFR-sparsity and DFR-TV schemes, we save on the complexity corresponding to solving the I2-I1 
and /2-TV optimization problems respectively. It is however hard to precisely give the order of complexity of solving the /2-^i 
and /2-TV optimization problems, as it is highly depend on the type of solvers. For some of the popular solvers like GPSR 
||33]| and NESTA ||35]| , the order of complexity is given as TO{NlogN), where T is the number of iterations and N is the 
resolution of the image [36|. Therefore, comparing to the DFR schemes we save a complexity of 2TO{NlogN). 

The complexity of the proposed scheme can be further reduced when a disparity value is estimated per block instead of per 
pixel. Fig. [7] compares the performance of our scheme with respect to the DFR schemes when disparity is estimated using per 
block with block size 4x4. Comparing Fig. [7] and Fig.^a) we see that the relative performances between the schemes remains 
approximately the same when the disparity image is estimated per block or per pixel. This confirms that the proposed scheme 
can easily adapt the granularity of disparity estimation without big penalty in order to meet the complexity requirements at 
the decoder. 



C Motion estimation performance 

We now illustrate the performance of our correlation estimation algorithm in video sequences. We first estimate motion 
vectors per blocks of 4 x 4 pixels with different sensing matrices <l>i 7^ <l>2 for each image. These vectors are then used to 
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Fig. 8. Illustration of the accuracy of motion field estimation in terms of image prediction quality for the Yosemite and Grove datasets. The quality of the 
predicted image 1 2 is compared with respect to 1 2 and /i . The prediction is carried out using the motion field estimated with block of pixels 4x4. 
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Fig. 9. Comparison of the quality of the predicted image I2 between the proposed, DFR-sparsity and DFR-TV schemes: (a) Yosemite dataset; (b) Grove 
dataset. The image prediction is carried out using the motion field that is estimated using block of pixels 4x4. 

predict the second image from the first image. Fig. [8] compares the predicted image I2 with the original images I2 and /i for 
Yosemite and Grove datasets respectively. It is clear that for a given measurement rate the predicted image I2 is closer to I2 
than /i which indicates that the motion between the images is efficiently captured by our correlation estimation algorithm. 
Similar experimental results are observed in the Mequon dataset. 

We then highlight the benefit of sampling the images with different sets of measurement matrices in Fig. [9j From Fig. [9] 
we see that the quality of the predicted image I2 is better when the images are sampled with different measurement matrices, 
compared to the case where the same sampling matrix is used for all images. This confirms the results shown for the disparity 
estimation performance. We finally compare our results to the DFR-sparsity and DFR-TV schemes that build the correlation 
model from (independently) reconstructed images based on minimizing the sparsity and TV priori respectively. Fig. [9ta) and 
Fig. Ob) show the comparison for the Yosemite and Grove datasets respectively. From Fig. [9] we see that for both datasets 
the proposed scheme performs better than the DFR-sparsity scheme at low rates and competes with the DFR-sparsity scheme 
at high rates, as observed in the disparity estimation study. Similar experimental findings are observed in the Mequon dataset. 
Furthermore, we see that the proposed scheme competes with the performance of DFR-TV scheme at low rate for the Grove 
dataset (see Fig. [3b)); as the Grove scene is textured with limited low frequency components, the TV prior in the reconstruction 
scheme results in poor reconstruction quality in the textured areas. Overall, the proposed scheme provides effective motion 
estimation results, while avoiding an order of computational complexity 2TO{NlogN) involved in the image reconstruction 
steps with the DFR-TV and DFR-sparsity schemes. 



D. Measurements quantization 

We briefly study here the performance of the correlation estimation algorithm when the measurements are affected by noise 
and in particular quantization noise. For the sake of simplicity we quantize the measurements using a uniform quantizer and 
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Fig. 10. Effect of quantization on the quality of the predicted image I2 in the Venus dataset. When the measurements are quantized we used BPDNp 071 
to reconstruct the images in DFR-TV based scheme. The image prediction is carried out using the dense disparity image estimated by solving Eq. f25> . 

we denote the quantized measurements as Yi and I2. We then estimate the correlation model by minimizing the energy in 
Eq. ([25l) using the quantized measurements Yi and I2. 

Fig. [TOl shows the effect of quantization on the disparity estimation performance when measurements are uniformly quantized 
using 2-, 3- and 4-bits. We use the disparity map to predict the second image in the Venus dataset. Interestingly we see that 
the 4-bit quantizer does not significantly affect the quality of the disparity image, as the degradation hardly reaches 0.5 dB in 
the quality of I2 at low to medium rates. As expected however the quality of the predicted image I2 decreases with increasing 
quantizer coarseness level for a fixed measurement rate. We then compare our results with the DFR-TV scheme that reconstructs 
the images by solving an optimization problem based on BPDNp in order to efficiently handle the quantization noise iJTll and 
then use the reconstructed images for disparity estimation. From Fig.[TOlwe see that the performance gap between the DFR-TV 
and compressed domain estimation is approximately the same in both the unquantized and quantized (i.e., 2-bit) scenarios. 
However, it should be noted that the DFR-TV scheme considers the nonlinearities due to quantization while reconstructing 
the images. Such effects are not considered in the proposed scheme. The solution of our scheme could also be improved by 
considering the quantization non-linearities but this problem is beyond the scope of this paper. 

We also study the effect of measurement quantization on the quality of the motion field. Fig. Ela) and Fig. fTTTb) shows the 
quality of the predicted image where the prediction is performed with motion vectors estimated from quantized measurements. 
The quality of the predicted image or equivalently the accuracy of motion estimation is reduced when the measurements 
are quantized as expected. Similarly to the case of disparity estimation, the influence of quantization is negligible when the 
measurements are quantized with a 4-bit quantizer. We also compare our results to the DFR-TV scheme that reconstructs 
the images by solving an optimization problem based on BPDNp and estimate motion from the reconstructed images. 
The performance of the DFR-TV scheme for the 2 -bit quantization scenario is shown in Fig. [TTl Interestingly, when the 
measurements are quantized we see that the proposed scheme competes with the DFR-TV scheme, because of the poor image 
reconstruction performance in the DFR-TV scheme when the measurements are coarsely quantized (i.e., 2-bit quantizer). 



E. Importance of correlation in joint reconstruction 

We finally propose to study the importance of accurate correlation estimation in a novel joint reconstruction algorithm (see 
Fig. [T]). We propose to reconstruct a pair of images /i and I2 by enforcing consistency with the compressed information and 
also with the estimated correlation model. A pair of image /i and I2 is reconstructed as a solution of the following constrained 
optimization problem: 



(^1,^2) =arg min (||V^*/i||i 

/l,/2 



IIV^*/2||i) s.t. ||yi-^i/i||2 = o,||F2- ^2/2112 = 0,11/2 



Ah\\l<e, 



(39) 



where ?/; is a redundant dictionary or an orthonormal basis in which the image is assumed to be sparse and ip* is the conjugate 
transpose of ip. From Eq. ([39l) it is clear that the images can be reconstructed independently if we solve the optimization 
without the last constraint Wh — ^^illi ^ ^- This corresponds to the independent reconstruction of sparse images in ip which 
agrees with the measurement information. By adding the last constraint, we impose that the pair of images also fit with the 
correlation model, in addition to the sparsity and data fidelity constraints. As a result the reconstruction quality for a given 
measurement rate is better when the images are reconstructed jointly than independently. The optimization problem for joint 
reconstruction can be re- written with proximity operators and solved efficiently using the parallel proximal algorithm (PPXA) 
proposed by Combettes et al 1381 (see 1391 for a similar solution in an asymmetric joint decoding scheme). 

We analyze the performance of our joint reconstruction scheme with a constraint imposed by the correlation estimated from 
linear measurements. In particular, we perform joint reconstruction experiments using the correlation model that is estimated 
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Fig. 11. Effect of using 4-, 3- and 2-bits quantizers on the image prediction quality for (a) Yosemite dataset and (b) Mequon dataset. When the measurements 
are quantized we used BPDNp IT/l to reconstruct the images in DFR-TV scheme. The image prediction is carried out using the motion field estimated by 
solving Eq. J25) 



32 
30 





^ 


<^ 


^ 


y 


-"" 


-" 




- 




— B— Joint reco: DFR-TV 
^ Joint reco: DFR-sparsity 

Joint reco: Proposed 

Indep reco 

Joint reco: originai 


- 











^ 


^ 


/ : 


— B— Joint reco: DFR-TV 

— -k — Joint reco: DFR-sparsity 

Joint reco: Proposed 

Indep reco 

Joint reco: originai 


- 









0.2 0.25 0.3 0.35 
Measurement rate per image 



0.2 0.25 0.3 0.35 

Measurement rate per image 



(a) 



(b) 



Fig. 12. Influence of disparity accuracy on the joint reconstruction performance for (a) Tsukuba and (b) Venus datasets. The joint reconstruction is carried 
out using the dense disparity image in the proposed and DFR-based schemes. The joint reconstruction performance is also compared to the independent 
reconstruction scheme. 



using a different measurement matrix for each image, i.e., <l>i ^ $2- We assume that the image is sparse in an orthonormal 
basis constructed using a wavelet transform. In our experiments, the parameter e in the optimization problem is selected based 
on a trial and error procedure that maximizes the reconstruction image quality /i and I2, and we set e = 14. We first compare 
our results to an independent reconstruction scheme that does not exploit the correlation between the images (i.e., the constraint 
1 1/2 — ^^i||2 ^ e is removed in Eq. (l39l)). Fig. fT2l compares the average reconstruction quality between the joint (denoted as 
Joint: Proposed) and independent reconstruction schemes for Tsukuba and Venus datasets. The joint reconstruction improves 
the reconstruction quality by 2 dB at low measurement rates and about 1 dB at high rates. We also observe in our experiments 
that the PSNR quality of the reconstructed images h and I2 is similar at a given measurement rate. The disparity estimate 
thus proves to be useful in improving the quality in the image reconstruction process. 

We then jointly reconstruct the images using a disparity estimated with DFR-sparsity and DFR-TV schemes. Fig. fT2l compares 
the quality of reconstructed image between the proposed and DFR schemes. We see that the disparity estimated from compressed 
measurements leads to a competitive performance with the disparity estimated by the DFR-sparsity scheme in terms of joint 
reconstruction quality. This is particularly obvious at low rate 0.05 where the DFR-sparsity scheme fails to accurately estimate 
the disparity. However, the quality of the reconstructed images is marginally penalized (i.e., 0.2 dB and 0.4 dB for the Tsukuba 
and Venus datasets respectively) compared to the reconstruction achieved when the disparity estimated by the DFR-TV scheme. 
Finally, we carry out the same experiments in a scenario where the images are jointly reconstructed using a correlation model 
that is estimated from the original images. This scheme serves as a benchmark for the joint reconstruction since the correlation 
is accurately known at the decoder. The corresponding results are denoted as Joint reco:original in Fig. [121 We see that the 
reconstruction quality achieved with the correlation estimated from compressed measurements converges to the performance 



benchmark when the measurement rate increases which further confirms the quahty of the disparity estimation. Finally, it 
should be noted that similar tendencies have been observed in joint reconstruction of video frames, but the corresponding 
results are omitted here due to space constraints. 

VI. Conclusions 

In this paper we have presented a framework for estimating the correlation between images given in the form of linear 
measurements without implementing explicit image reconstruction steps. We have proposed a linear representation of disparity 
and motion models and show that the correlation can be estimated in the compressed domain, thanks to the distance preserving 
property of the sensing matrix. The correlation is estimated by solving a regularized energy model that enforces consistency with 
the measurements and smoothness of the correlation information. Extensive experimental results demonstrate that the proposed 
methodology provides a good estimation of disparity or motion fields in different natural and synthetic image datasets, especially 
when images are sampled with different measurement matrices. We also show that our correlation estimation solution competes 
with the correlation estimation solutions that reconstructs images a priori, but becomes clearly advantageous due to its lower 
computational complexity. The correlation estimation from compressed measurements thus provides an effective solution for 
distributed scene analysis or coding applications in low complexity sensor networks. 
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